% Part 1: Bayesian Theory Foundation
% This file is included by bayesian_vs_3dvar_cn_main.tex

\section{大气反演的贝叶斯理论基础}

\subsection{贝叶斯定理在大气反演中的应用}
根据贝叶斯定理,在给定观测 $\mathbf{y}^{\text{obs}}$ 条件下,通量状态 $\mathbf{x}$ 的后验概率分布为
\begin{equation}
p(\mathbf{x} | \mathbf{y}^{\text{obs}}) = \frac{p(\mathbf{y}^{\text{obs}} | \mathbf{x}) p(\mathbf{x})}{p(\mathbf{y}^{\text{obs}})},
\end{equation}
其中 $p(\mathbf{y}^{\text{obs}} | \mathbf{x})$ 为似然函数,$p(\mathbf{x})$ 为先验分布,$p(\mathbf{y}^{\text{obs}})$ 为观测的边缘概率(作为归一化常数)。

{\color{red}\textbf{[贝叶斯定理核心]：后验概率 = 似然度 × 先验概率 / 归一化常数。这个公式告诉我们如何用观测数据更新对参数的信念。似然度衡量模型与数据的拟合程度，先验体现已有知识，两者的乘积经归一化后给出更新后的后验分布。}}

\subsection{前向模型}
大气观测与地表通量的关系通过线性输送算子 $\mathbf{H}$ 联系:
\begin{equation}
\mathbf{y}^{\text{obs}} = \mathbf{H}\mathbf{x} + \boldsymbol{\varepsilon},
\end{equation}
其中 $\mathbf{y}^{\text{obs}} \in \mathbb{R}^m$ 为观测向量,$\mathbf{x} \in \mathbb{R}^n$ 为通量状态向量,$\mathbf{H} \in \mathbb{R}^{m \times n}$ 为源-受体关系矩阵,$\boldsymbol{\varepsilon}$ 为观测误差。

{\color{red}\textbf{[输送算子的物理意义]：矩阵 $\mathbf{H}$ 的每个元素 $H_{ij}$ 表示第 $j$ 个网格的单位通量对第 $i$ 个观测点混合比的贡献。在 FLEXPART 中，这通过后向拉格朗日粒子的驻留时间计算得到，精确捕捉大气输送过程。}}

对于非全球反演域,需区分嵌套域内外的通量贡献:
\begin{equation}
\mathbf{y}^{\text{obs}} = \mathbf{H}^{\text{nest}}\mathbf{x}^{\text{nest}} + \mathbf{H}^{\text{out}}\mathbf{x}^{\text{out}} + \mathbf{H}^{\text{bg}}\mathbf{y}^{\text{bg}} + \boldsymbol{\varepsilon},
\end{equation}
其中 $\mathbf{H}^{\text{nest}}$ 为嵌套域输送算子,$\mathbf{H}^{\text{out}}$ 为域外贡献,$\mathbf{H}^{\text{bg}}\mathbf{y}^{\text{bg}}$ 为初始混合比贡献。

\subsection{高斯假设下的代价函数}

\subsubsection{从贝叶斯定理到代价函数的推导}
本小节展示第2.1节贝叶斯定理与代价函数的内在联系。在高斯假设下,贝叶斯后验分布可转化为等价的优化问题。

假设先验和观测误差服从高斯分布:
\begin{align}
p(\mathbf{x}) &= \mathcal{N}(\mathbf{x}_b, \mathbf{B}) = \frac{1}{(2\pi)^{n/2}|\mathbf{B}|^{1/2}} \exp\left[-\frac{1}{2}(\mathbf{x} - \mathbf{x}_b)^T\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}_b)\right], \\
p(\mathbf{y}^{\text{obs}} | \mathbf{x}) &= \mathcal{N}(\mathbf{H}\mathbf{x}, \mathbf{R}) = \frac{1}{(2\pi)^{m/2}|\mathbf{R}|^{1/2}} \exp\left[-\frac{1}{2}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}})^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}})\right].
\end{align}

根据第2.1节的贝叶斯定理,后验分布为
\begin{equation}
p(\mathbf{x} | \mathbf{y}^{\text{obs}}) \propto p(\mathbf{y}^{\text{obs}} | \mathbf{x}) p(\mathbf{x}).
\end{equation}

将高斯分布代入并取负对数,得到
\begin{align}
-\ln p(\mathbf{x} | \mathbf{y}^{\text{obs}}) &= -\ln p(\mathbf{y}^{\text{obs}} | \mathbf{x}) - \ln p(\mathbf{x}) + \text{const} \\
&= \frac{1}{2}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}})^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}}) \nonumber \\
&\quad + \frac{1}{2}(\mathbf{x} - \mathbf{x}_b)^T\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}_b) + \text{const}.
\end{align}

{\color{red}\textbf{[对数变换的关键作用]：取负对数将概率密度的\textbf{乘积}转化为\textbf{和}，将\textbf{最大化后验概率}(MAP) 转化为\textbf{最小化代价函数}。高斯分布的负对数恰好是二次型，使问题变为凸优化。归一化常数 $(2\pi)^{n/2}|\mathbf{B}|^{1/2}$ 与 $\mathbf{x}$ 无关，在优化中可忽略。}}

因此,最大后验估计等价于最小化代价函数
\begin{equation}
J(\mathbf{x}) = \frac{1}{2}(\mathbf{x} - \mathbf{x}_b)^T\mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}_b) + \frac{1}{2}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}})^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}}),
\end{equation}
其中 $\mathbf{x}_b$ 为先验通量,$\mathbf{B}$ 为先验误差协方差矩阵,$\mathbf{R}$ 为观测误差协方差矩阵。第一项衡量偏离先验估计的代价,第二项衡量模拟与观测的不匹配程度。

{\color{red}\textbf{[代价函数的平衡]：这两项体现了贝叶斯方法的核心哲学：在先验约束与观测拟合之间寻求平衡。$\mathbf{B}^{-1}$ 和 $\mathbf{R}^{-1}$ 作为精度矩阵，自动调节两项的相对权重。当观测不确定性大($\mathbf{R}$ 大)时，解更接近先验；当先验不确定性大($\mathbf{B}$ 大)时，解更贴近观测。}}

\subsubsection{高斯假设的关键性}
上述推导严格依赖高斯假设。对于非高斯情况:
\begin{itemize}
\item \textbf{非高斯先验}(如对数正态、伽马分布): 代价函数不再是二次型,需使用非线性优化或MCMC
\item \textbf{非高斯似然}(如Student-t、混合高斯): $-\ln p(\mathbf{y}^{\text{obs}} | \mathbf{x})$ 不再简化为二次型
\item \textbf{多峰后验}: MAP估计可能陷入局部极值,无法代表完整后验分布
\end{itemize}

这些情况下,第2.1节的通用贝叶斯框架依然成立,但必须采用MCMC等采样方法探索完整后验分布,而非简单地最小化代价函数。

\subsection{解析解}
代价函数的梯度为
\begin{equation}
\nabla J(\mathbf{x}) = \mathbf{B}^{-1}(\mathbf{x} - \mathbf{x}_b) + \mathbf{H}^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{x} - \mathbf{y}^{\text{obs}}).
\end{equation}
令 $\nabla J(\mathbf{x}) = 0$,可得最优估计
\begin{equation}
\mathbf{x}_a = \mathbf{x}_b + \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T + \mathbf{R})^{-1}(\mathbf{y}^{\text{obs}} - \mathbf{H}\mathbf{x}_b),
\end{equation}
对应的后验误差协方差矩阵为
\begin{equation}
\mathbf{A} = (\mathbf{B}^{-1} + \mathbf{H}^T\mathbf{R}^{-1}\mathbf{H})^{-1} = \mathbf{B} - \mathbf{B}\mathbf{H}^T(\mathbf{H}\mathbf{B}\mathbf{H}^T + \mathbf{R})^{-1}\mathbf{H}\mathbf{B}.
\end{equation}

{\color{red}\textbf{[计算复杂度]：解析解需要求逆 $m \times m$ 矩阵($m$ 为观测数)，复杂度为 $\mathcal{O}(m^3)$。对于稀疏观测($m \ll n$)效率高，但对于密集观测或卫星数据($m \sim 10^5$-$10^6$)，直接求逆不可行，需采用共轭梯度等迭代方法。}}

\subsection{控制变量变换}
为改善数值稳定性,引入 $\chi$ 空间变换
\begin{equation}
\chi = \mathbf{B}^{-1/2}(\mathbf{x} - \mathbf{x}_b),
\end{equation}
使控制变量在统计上无相关且方差为 1。在 $\chi$ 空间中,代价函数简化为
\begin{equation}
J(\chi) = \frac{1}{2}\chi^T\chi + \frac{1}{2}(\mathbf{H}\mathbf{B}^{1/2}\chi - \mathbf{d})^T\mathbf{R}^{-1}(\mathbf{H}\mathbf{B}^{1/2}\chi - \mathbf{d}),
\end{equation}
其中 $\mathbf{d} = \mathbf{y}^{\text{obs}} - \mathbf{H}\mathbf{x}_b$ 为创新向量。

{\color{red}\textbf{[预条件化的作用]：$\chi$ 空间变换将先验协方差矩阵 $\mathbf{B}$ 的本征值谱归一化，显著改善问题的条件数。这类似于将一个扁平的椭球变换为圆球，使各方向的优化难度均等。对于共轭梯度法，条件数从 $\kappa(\mathbf{B}) \sim 10^6$ 降至 $\kappa(\mathbf{I}) = 1$，收敛速度提升数量级。}}

\subsection{矩阵构造概览}

FLEXINVERT 系统通过三个关键矩阵实现贝叶斯框架：

\textbf{先验协方差矩阵 $\mathbf{B}$}: 采用指数空间相关模型(海陆差异化相关长度)与时间相关模型,通过本征值分解实现高效存储($\mathcal{O}(KN)$ 而非 $\mathcal{O}(N^2)$)。保留主要本征模态 $K \ll N$,消除近零本征值改善条件数。

\textbf{观测误差协方差矩阵 $\mathbf{R}$}: 对角结构,误差包含仪器精度、代表性误差和输送模型不确定性。根据观测类型(自由对流层、近海岸)自适应调整误差权重。

\textbf{输送算子 $\mathbf{H}$}: 通过 FLEXPART 后向拉格朗日粒子追踪计算源-受体关系。$H_{ij}$ 表示第 $j$ 个网格单位通量对第 $i$ 个观测点的贡献。采用变分辨率聚合策略,在高敏感度区域保持细分辨率。

数值求解使用解析解($N < 10^4$)、共轭梯度法(中等规模)或拟牛顿法 M1QN3($N \sim 10^6$)。

详细的矩阵构造方法、数学公式和实现细节请参见附录~\ref{app:flexinvert-matrices}(第~\pageref{app:flexinvert-matrices}页)。

